//----------------------------------*-C++-*----------------------------------//
/*!
 * \file   SlabtranInput.cc
 * \author Jeremy Roberts
 * \date   04/10/2011
 * \brief  Member definitions of class SlabtranInput
 * \note   Copyright (C) 2011 Jeremy Roberts
 */
//---------------------------------------------------------------------------//
// $Rev:: 141                                           $:Rev of last commit
// $Author:: j.alyn.roberts@gmail.com                   $:Author of last commit
// $Date:: 2011-09-14 19:28:23 +0000 (Wed, 14 Sep 2011) $:Date of last commit
//---------------------------------------------------------------------------//
#include <iostream>
#include <vector>

#include "tran/SlabtranInput.hh"
#include "../utilities/GenException.hh"

namespace slabtran
{

//---------------------------------------------------------------------------//
// CONSTRUCTOR
//---------------------------------------------------------------------------//
SlabtranInput::SlabtranInput(string S)
: InputXML(S)
{
	// control defaults
	d_ptype 	= -1;
	d_stype 	= 0;
	d_quadord 	= 4;
	d_epss 		= 0.123;
	d_epsk 		= 0.124;
	d_maxit 	= 999;
	d_angord    = 2;
	d_numel 	= 1;
	// output defaults
	d_printout    = -1;
	d_plotout     = -1;
	d_outfile     = "n/a";
	d_fluxfile    = "n/a";
	// material defaults
	numg        = -1;
	numm        = -1;
}

//---------------------------------------------------------------------------//
// DESTRUCTOR
//---------------------------------------------------------------------------//
SlabtranInput::~SlabtranInput()
{
	return;
}

////---------------------------------------------------------------------
//// Access element contents for a single element

// Check element index
void SlabtranInput::checkElIdx( int e )
{
	if ( e < 0 || e >= elements.size() )
	{
		throw GenException(__LINE__,__FILE__,
				"Element index out of bounds.");
	}
	return;
}

// Return my identification
int SlabtranInput::id( int e )
{ 
	checkElIdx(e);
	return elements[e].d_id; 
}

// Return my description
string SlabtranInput::desc( int e )
{ 
	checkElIdx(e);
	return elements[e].d_desc; 
}

// Return source
double SlabtranInput::src( int e, int i, int j, int k, int g )
{ 
	checkElIdx(e);
	return elements[e].d_src[i][j][k][g]; 
}

// Return coarse mesh count in in x/y/z
int SlabtranInput::ncm( int e, int dim )
{ 
	checkElIdx(e);
	return elements[e].d_ncm[dim]; 
}

// Return coarse mesh count in in x/y/z
double SlabtranInput::cm( int e, int xyz, int dim )
{ 
	checkElIdx(e);
	return elements[e].d_cm[xyz][dim]; 
}

// Return fine mesh count in x/y/z
int SlabtranInput::nfm( int e, int xyz, int dim )
{ 
	checkElIdx(e);
	return elements[e].d_nfm[xyz][dim]; 
}

// Return fine mesh count in x/y/z
int SlabtranInput::mt( int e, int i, int j, int k)
{
	checkElIdx(e);
	return elements[e].d_mt[i][j][k];
}

// Return SP to element
SlabtranInput::SP_element SlabtranInput::element( int e )
{
	checkElIdx(e);
	SlabtranElement *f = &elements[e];
	SP_element spf(f);
	return spf;
}

//---------------------------------------------------------------------------//
/*!
 * \brief This processes the global input.
 *
 */
bool SlabtranInput::processFile()
{

	//--------------------------------------------------------------------------
	// get problem name

	xmlNodePtr cur = getXpathElement("/slabtraninput");
	this->d_name = (getXpathContent(cur,"name"));
	cout << " inp.name = " << this->d_name << endl;

	//--------------------------------------------------------------------------
	// move on to control

	cur      = getXpathElement("/*/control");
	d_ptype = atoi( ( getXpathContent(cur,"ptype") ) );
	d_stype = atoi( ( getXpathContent(cur,"stype") ) );
	d_spacedim = atoi( ( getXpathContent(cur,"spacedim") ) );
	d_numel = atoi((getXpathContent(cur,"numel")));
	if ( nodeExists( cur, "quadord" ) and d_stype==1 )
		d_quadord = atoi(getXpathContent( cur, "quadord" ));
	else if ( d_stype== 1 )
		d_quadord = 4;
	else
		d_quadord = 1; // have to be careful about diffusion
	if ( nodeExists( cur, "epsk" ) )
		d_epsk = atof((getXpathContent(cur,"epsk")));
	if ( nodeExists( cur, "epss" ) )
		d_epss = atof((getXpathContent(cur,"epss")));
	if ( nodeExists( cur, "maxit" ) )
		d_maxit = atoi((getXpathContent(cur,"maxit")));
	if ( nodeExists( cur, "keff" ) )
		d_keff = atof((getXpathContent(cur,"keff")));
	if ( nodeExists( cur, "angord" ) )
		d_angord = atoi((getXpathContent(cur,"angord")));

	//--------------------------------------------------------------------------
	// move on to output

	cur      = getXpathElement("/*/output");
	d_printout = atoi((getXpathContent(cur,"printout")));
	d_plotout  = atof((getXpathContent(cur,"plotout")));
	if ( nodeExists( cur, "outfile" ) )
		d_outfile = getXpathContent(cur,"outfile");
	if ( nodeExists( cur, "fluxfile" ) )
		d_fluxfile = getXpathContent(cur,"fluxfile");

	//--------------------------------------------------------------------------
	// materials  --- eventually, this stuff should be stored in a mat class

	cur  = getXpathElement("/*/materials");
	info = getXpathContent(cur,"info");
	numm = atoi( ( getXpathContent(cur,"numm") ) );
	numg = atoi( ( getXpathContent(cur,"numg") ) );

	matname.resize(numm);
	if ( d_stype == 0 ) dc.resize(numm);
	if ( d_stype == 1 ) st.resize(numm);
	sr.resize(numm);
	ab.resize(numm);
	ns.resize(numm);
	chi.resize(numm);
	ss.resize(numm);
	for(int i = 0; i < numm; ++i)
	{
		if ( d_stype == 0 ) dc[i].resize(numg);
		if ( d_stype == 1 ) st[i].resize(numg);
		sr[i].resize(numg);
		ab[i].resize(numg);
		ns[i].resize(numg);
		chi[i].resize(numg);
		ss[i].resize(numg);
		for(int j = 0; j < numg; ++j)
			ss[i][j].resize(numg);
	}

	xmlNodeSetPtr mats = getXpathElements("/*/materials/mat");
	xmlNodeSetPtr grps;
	xmlNodeSetPtr sigs;

	cout << "number of materials " << mats->nodeNr << endl;
	if ( mats->nodeNr != numm )
	{
		cout << "Input Error: <mat> count doesn't match numm" << endl;
		return false;
	}

	for ( int m = 0; m < numm; m++ )
	{

		cout << " ---------------------------- " << endl;
		cout << " mat number = " << m << endl;
		cout << " mat name   = "
				<< getXpathContent(mats->nodeTab[m],"matname") << endl;
		matname[m] = getXpathContent(mats->nodeTab[m],"matname");
		cout << " node name = " << mats->nodeTab[m]->name << endl;
		cur = mats->nodeTab[m];
		grps = getXpathElements( cur, "group" );

		if ( grps->nodeNr != numg )
		{
			cout << "Input Error: <mat> number " << m
					<< " group count doesn't match numg" << endl;
			return false;
		}

		for ( int g = 0; g < numg; g++ )
		{
			cout << "       ---------------------------- " << endl;
			cout << "       group number = " << g << endl;
			cout << "       node name = " << grps->nodeTab[g]->name << endl;
			cur =  grps->nodeTab[g];
			if ( d_stype == 0 ) dc[m][g] = atof( getXpathContent( cur, "D" ) );
			if ( d_stype == 1 ) st[m][g] = atof( getXpathContent( cur, "sT" ) );
			sr[m][g] = atof( getXpathContent( cur, "sR" ) );
			ns[m][g] = atof( getXpathContent( cur, "nsF" ) );
			chi[m][g] = atof( getXpathContent( cur, "chi" ) );
			sigs = getXpathElements( cur, "sS" );
			if ( sigs->nodeNr != numg )
			{
				cout << "Input Error: <mat> number " << m << " <group> number "
						<< g << " sS count doesn't match numg" << endl;
				return false;
			}
			for ( int gg = 0; gg < numg; gg++ )
			{
				cur = sigs->nodeTab[gg];
				// ss[material][to group][from group]
				ss[m][g][gg] = atof( (const char*)cur->xmlChildrenNode->content );
			}   // end gg loop
		}       // end g loop

		// compute absorption cross-section
		for ( int g = 0; g < numg; g++ )
		{
			scalar sum = 0;
			for( int gg = 0; gg < numg; ++gg )
			{
				if ( gg != g ) sum = sum + ss[m][gg][g];
				// othewise, we are ignoring inscatter
				ab[m][g] = sr[m][g] - sum;
			} // end gg loop 
		}     // end g loop

	}         // end m loop

	//--------------------------------------------------------------------------
	// geometry & elements
	xmlNodeSetPtr els = getXpathElements("/*/geometry");
	if ( els->nodeNr != d_numel )
	{
		cout << "Input Error: <geometry> count " << els->nodeNr
				<< " doesn't match numel " << d_numel << endl;
		return false;
	}

	elements.resize(d_numel);

	xmlNodeSetPtr cm;

	for ( int e = 0; e < d_numel; e++ )
	{

		cout << " ---------------------------- " << endl;
		cout << " element number      = " << e << endl;
		cout << " element discription = "
				<< getXpathContent(els->nodeTab[e],"description") << endl;

		cur = els->nodeTab[e];
		elements[e].N = d_spacedim;
		elements[e].d_id = e;
		elements[e].d_desc = getXpathContent(els->nodeTab[e],"description");

		elements[e].d_bc.resize(2*d_spacedim);
		elements[e].d_bc[LEFT]  = atof( getXpathContent( cur, "bcl" ) );
		elements[e].d_bc[RIGHT] = atof( getXpathContent( cur, "bcr" ) );
		if ( d_spacedim > 1 )
		{
			elements[e].d_bc[BOTTOM]  = atof( getXpathContent( cur, "bcb" ) );
			elements[e].d_bc[TOP] = atof( getXpathContent( cur, "bct" ) );
		}
		if ( d_spacedim > 2 )
		{
			elements[e].d_bc[NORTH]  = atof( getXpathContent( cur, "bcn" ) );
			elements[e].d_bc[SOUTH] = atof( getXpathContent( cur, "bcs" ) );
		}

		elements[e].d_ncm.resize(3);
		elements[e].d_ncm[I] = atof( getXpathContent( cur, "nxcm" ) );
		if ( d_spacedim > 1 )
			elements[e].d_ncm[J] = atof( getXpathContent( cur, "nycm" ) );
		else
			elements[e].d_ncm[J] = 1;
		if ( d_spacedim > 2 )
			elements[e].d_ncm[K] = atof( getXpathContent( cur, "nzcm" ) );
		else
			elements[e].d_ncm[K] = 1;
		elements[e].d_cm.resize(3); // resize first dimension of cm edges
		elements[e].d_nfm.resize(3); // resize first dimension of fm counts

		// DIM 1 ----------------- perhaps functionalize this w/ for a loop i=1,spacedim

		// xcm
		cur = getXpathElement( els->nodeTab[e], "xcm" );
		cm = getXpathElements( cur, "e" );
		if ( cm->nodeNr != elements[e].d_ncm[I] + 1 )
		{
			std::cout << " issue?? " << std::endl;
		}
		elements[e].d_cm[I].resize(elements[e].d_ncm[I]+1);
		elements[e].d_cm[J].resize(1);  // warning: hard wired 1d for now
		elements[e].d_cm[K].resize(1);
		for ( int i = 0; i < elements[e].d_ncm[I]+1; i++ )
		{
			cur = cm->nodeTab[i];
			elements[e].d_cm[I][i] = atof( (const char*)cur->xmlChildrenNode->content );
		}
		// nxfm
		cur = getXpathElement( els->nodeTab[e], "nxfm" );
		cm = getXpathElements( cur, "n" );
		if ( cm->nodeNr != elements[e].d_ncm[I] )
		{
			cout << "Input Error: <element> number " << e << " nfxm entry count "
					<< cm->nodeNr << " doesn't match nxcm = "
					<< elements[e].d_ncm[I] << endl;
			return false;
		}
		elements[e].d_nfm[I].resize(elements[e].d_ncm[I]);
		elements[e].d_nfm[J].resize(elements[e].d_ncm[J]); // warning: hard wired 1d for now
		elements[e].d_nfm[K].resize(elements[e].d_ncm[K]);
		for ( int i = 0; i < elements[e].d_ncm[I]; i++ )
		{
			cur = cm->nodeTab[i];
			elements[e].d_nfm[I][i] = atoi( (const char*)cur->xmlChildrenNode->content );
		}
		elements[e].d_nfm[J][0] = 1;
	    elements[e].d_nfm[K][0] = 1;
		// mt

		// allocate. it might be better to have classes for these matrix quantities!
		elements[e].d_mt.resize( elements[e].d_ncm[I] );
		for ( int i = 0; i < elements[e].d_ncm[I]; i++ )
		{
			elements[e].d_mt[i].resize( elements[e].d_ncm[J] );
			for ( int j = 0; j < elements[e].d_ncm[J]; j++ )
			{
				elements[e].d_mt[i][j].resize( elements[e].d_ncm[K] );
			}
		}

		cur = getXpathElement( els->nodeTab[e], "matid" );
		cm = getXpathElements( cur, "m" );
		// warning: must account for multi-d
//		if ( cm->nodeNr != elements[e].d_nxcm  )
//		{
//			cout << "Input Error: <element> number " << e << " matid entry count "
//					<< cm->nodeNr << " doesn't match nxcm = "
//					<< elements[e].d_nxcm << endl;
//			return false;
//		}
		// warning: hardcoding 1-d read for now.
		for ( int i = 0; i < elements[e].d_ncm[I]; i++ )
		{
			cur = cm->nodeTab[i];
			elements[e].d_mt[i][0][0] = atoi( (const char*)cur->xmlChildrenNode->content );
		}

		// source

		// allocate. it might be better to have classes for these matrix quantities!
		elements[e].d_src.resize( elements[e].d_ncm[I] );
		for ( int i = 0; i < elements[e].d_ncm[I]; i++ )
		{
			elements[e].d_src[i].resize( elements[e].d_ncm[J] );
			for ( int j = 0; j < elements[e].d_ncm[J]; j++ )
			{
				elements[e].d_src[i][j].resize( elements[e].d_ncm[K] );
				for ( int k = 0; k < elements[e].d_ncm[J]; k++ )
				{
					elements[e].d_src[i][j][k].resize( numg );
				}
			}
		}

		// we're still in the geometry loop; get this guy's source
		cur  = getXpathElement( els->nodeTab[e], "source" );
		xmlNodeSetPtr sgrps = getXpathElements( cur, "sgroup" );
		cout << " num sgroups = " << sgrps->nodeNr << endl;
		if ( sgrps->nodeNr != numg )
		{
			std::cout << " number sgroups != numg " << std::endl;
		}
		for ( int g = 0; g < numg; g++ )
		{
			cur = sgrps->nodeTab[g];
			cm  = getXpathElements(cur,"s");
			if ( cm->nodeNr != elements[e].d_ncm[I] )
			{
				std::cout << " source has wrong number of s elements = " << cm->nodeNr << std::endl;
 			}
			for ( int i = 0; i < elements[e].d_ncm[I]; i++ )
			{
				cur = cm->nodeTab[i];
				elements[e].d_src[i][0][0][g] = atof( (const char*)cur->xmlChildrenNode->content );
			}
		}

	}

	return true;
}

//---------------------------------------------------------------------------//
/*!
 * \brief This echos the input.
 *
 */
void SlabtranInput::echoInput()
{
	cout << " echo of user input... " << endl;
	cout << endl;
	cout << " Problem name:  " << name() << endl;
	cout << " --------------------------------------------- " << endl;
	cout << " Control" << endl;
	cout << " --------------------------------------------- " << endl;
	cout << "          ptype: " << ptype() << endl;
	cout << "          stype: " << stype() << endl;
	cout << "        quadord: " << quadord() << endl;
	cout << "           epss: " << epss() << endl;
	cout << "           epsk: " << epsk() << endl;
	cout << "          maxit: " << maxit() << endl;
	cout << "         angord: " << angord() << endl;
	cout << "          numel: " << numel() << endl;
	cout << " --------------------------------------------- " << endl;
	cout << " Output" << endl;
	cout << " --------------------------------------------- " << endl;
	cout << "       printout: " << d_printout << endl;
	cout << "        plotout: " << d_plotout << endl;
	cout << "        outfile: " << d_outfile << endl;
	cout << "       fluxfile: " << d_fluxfile << endl;
	cout << " --------------------------------------------- " << endl;
	cout << " Materials" << endl;
	cout << " --------------------------------------------- " << endl;
	cout << "           info: " << info << endl;
	cout << "           numm: " << numm << endl;
	cout << "           numg: " << numg << endl;
	for ( int m = 0; m < numm; m++ )
	{
		cout << "       material: " << m << endl;
		cout << "        matname: " << matname[m] << endl;
		for ( int g = 0; g < numg; g++ )
		{
			cout << "          group: " << g << endl;
			cout << "                   ";
			if ( stype() == 0 )
			{
				printf (" %10.6f ", dc[m][g] );
			}
			if ( stype() == 1 )
			{
				printf (" %10.6f ", st[m][g] );
			}
			printf (" %10.6f ", sr[m][g] );
			printf (" %10.6f ", ab[m][g] );
			printf (" %10.6f ", ns[m][g] );
			printf (" %10.6f ", chi[m][g] );
			for ( int gg = 0; gg < numg; gg++ )
			{
				//                ss[material][from][to]
				printf (" %10.6f ", ss[m][g][gg] );
			}
			cout << endl;

		}
		cout << endl;
	}
	cout << " --------------------------------------------- " << endl;
	cout << " ELEMENTS" << endl;
	cout << " --------------------------------------------- " << endl;
	for ( int e = 0; e < numel(); e++ )
	{
		cout << endl;
		cout << "        element: " << e << endl;
		cout << "    description: " << elements[e].desc()      << endl;
		cout << "            bcl: " << elements[e].d_bc[LEFT]  << endl;
		cout << "            bcr: " << elements[e].d_bc[RIGHT] << endl;
		cout << "           nxcm: " << elements[e].d_ncm[I]      << endl;
		cout << "            xcm: ";
		for ( int i = 0; i < elements[e].d_ncm[I]+1; i++ )
			printf (" %8.4f ", elements[e].d_cm[I][i] );
		cout << endl;
		cout << "           nxfm:   ";
		for ( int i = 0; i < elements[e].d_ncm[I]; i++ )
			printf (" %8i ", elements[e].d_nfm[I][i] );
		cout << endl;
		cout << "             mt:   ";
		for ( int i = 0; i < elements[e].d_ncm[I]; i++ )
			printf (" %8i ", elements[e].mt(i) );
		cout << endl;
		cout << "            src:   " << endl;
		for ( int g = 0; g < numg; g++ )
		{
			cout << "          group: " << g << endl;
			cout << "                   ";
			for ( int i = 0; i < elements[e].d_ncm[I]; i++ )
				printf ("%8.4f ", elements[e].d_src[i][0][0][g] );
			cout << endl;
		}
	}

	return;
}

//---------------------------------------------------------------------------//
/*!
 * \brief This gives a code announcement.
 *
 */
void SlabtranInput::sayHello()
{
	//      _          _      _
	// ___ | |  __ _  | |__  | |_   _ _   __ _   _ _
	//(_-< | | / _` | | '_ \ |  _| | '_| / _` | | ' \
	///__/ |_| \__,_| |_.__/  \__| |_|   \__,_| |_||_|
	cout << "|-------------------------------------------------|" << endl;
	cout << "|                 S L A B T R A N                 |" << endl;
	cout << "|  a simple 1-d transport code for slab geometry  |" << endl;
	cout << "|                                                 |" << endl;
	cout << "|  Version 0.1, 04/10/2011                        |" << endl;
	cout << "|-------------------------------------------------|" << endl;
}

} // end namespace slabtran

//---------------------------------------------------------------------------//
//                 end of SlabtranInput.cc
//---------------------------------------------------------------------------//

